Growth-Collapse Cycles of a Bose-Einstein Condensate with Attractive Interactions 
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A Bose-Einstein condensate of atoms with attractive interactions exhibits growth and collapse 
cycles, when it is fed by a thermal cloud. Recently this phenomenon has been directly observed in a 
trapped 7 Li gas. We offer a quantitative explanation of the data, on the basis of a model proposed 
earlier, ft is shown that the condensate wave function acquires a chaotic component after the first 
■ collapse, indicating superfluid turbulence. 
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I. INTRODUCTION 

o 

In an elegant experiment jLl, Hulet and his team at Rice University have observed the growth and collapse of a 
Bose-Einstein condensate of magnetically trapped 7 Li atoms, whose interactions are predominantly attractive. This 
fascinating system had been created and studied || by the same group some time ago; but the instability of the 
system has not been directly observed until now. The purpose of this paper is to give a quantitative explanation of 
the data, based on a theoretical model proposed earlier Q ||. We refer to these references for earlier literature. 

Because of the attractive interactions, 7 Li atoms in free space would solidify at low temperatures. When confined 
in an external potential, however, the atoms are maintained in a gaseous state by their zero-point motion, as long as 
their number is small. Above a critical number, however, the attractive interactions overcome the kinetic energy, and 
the system collapses into a state of high density. This instability is similar to that of a white dwarf star whose mass is 
greater than the Chandrasekhar limit, in which the gravitational attraction overwhelms the zero-point pressure of the 
, Fermi gas of electrons in the star. Unlike in a white dwarf star, however, the particle number here is not conserved, for 
the atoms can be expelled from the magnetic trap due to two- and three-body collisions, which become increasingly 
important at higher densities. During the collapse, the particle number suddenly drops to a subcritical value, only 
to grow back subsequently, if the condensate is fed by a surrounding thermal cloud. Thus, the condensate undergoes 
cycles of growth and collapse, and in this respect the system is more akin to an electric buzzer than a white dwarf 
star. 

Let us first ignore the gain and loss mechanisms. To model the system, we treat the condensate in the mean-field 
approximation, an assume that the condensate wave function ip satisfies a Gross-Pitaevskii equation, or nonlinear 
C$ ' Schrodinger equation (NLSE): 
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U = ^, (1) 
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where a is a negative scattering length, and the external potential is taken to be harmonic: 



V(r) = imwV. (2) 
This defines a characteristic length do, the width of the unperturbed ground-state wave function: 



d = d—. (3) 
V mu> 

The number of condensate particles enters through the normalization 

N= I d 3 r\4>\ 2 . (4) 
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a = —1.45 nm, 
do ~ 3.16 /Ltm, 

luk 908 s- 1 . (6) 

Actually, the magnetic trap was not quite spherically symmetric, and u) here represents the geometric mean of the 
frequencies. 

In the absence of an external potential in spatial dimensions D > 2, the NLSE with attractive interactions (Uo > 0) 
exhibits "self- focusing" |fj), in which the system collapses to a state of locally infinite density in a finite time. The 
cause of this instability is the term — i[/o|V'| 2 in the Hamiltonian ||, which has no lower bound as the density \ip\ 2 
increases. In the presence of an external trap V(r), however, there exists an energy barrier against collapse. As a 
function of l^l 2 , the form of this barrier depends on the location r, and its height is smallest at r = 0, where the 
external potential is weakest. As predicted in Q, and verified through numerical calculations in the energy barrier 
at r = is breached when N > N c , with 

N c = 0.574-^ = 1250, (7) 
\a\ 

and this triggers the collapse of the system. At the center of the trap, there appears a black hole into which particles 
are drawn, and they will be drained from the system when loss mechanisms are taken into account. Of course, particles 
can surmount the energy barrier through quantum tunneling; even when N < N c . but this is a slow process that can 
be neglected compared to loss through collisions. 

What happens to the particles that go through the energy barrier? According to the Hamiltonian [5] their density 
would continue to increase indefinitely; but, of course, once the density begins to grow, effects so far neglected 
will become important. An example is solidification, which might be simulated by adding high-order terms in the 
Hamiltonian, such as c|?/>| 6 , with c > 0. Such a term in the NLSE will lead to the formation of droplets in a first-order 
phase transition 0. However, such effects are overshadowed by loss though collisions, and will be ignored. 



The corresponding Hamiltonian is given by 



H 



d A r 



In the Rice experiment, the parameters have the values 



II. GAIN AND LOSS MECHANISMS 



Ketterlc and his team at MIT has measured the rate of growth of a 23 Na condensate fed by a thermal cloud, and 
fitted the results to the following formula H : 



dN 



N 



(8) 



with Co = (50 ms) _1 , N cq w5x 10 5 . The factor in square brackes represents a saturation effect. We shall ignore it by 
assuming N << N eq for our case. The growth rate Co should be proportional to the elastic scattering cross section, 
and hence to the square of the scattering length. Since the scattering lengths for 23 Na and 7 Li are 2.75 nm and —1.45 
nm respectively, we take for 7 Li 



co 



1.45 

2J5 



(50 ms)- 1 . 



(9) 



The loss rate ^(Af) due to two and three-body collisions have been calculated by Dodd et al J|/ ; who report their 
results in the form 



R(N)=A / d 3 r\i[>\ 4 + B / d 3 r|</f, 



A = 1.2 x 10~ 14 cm 3 s -1 , 
B = 2.6 x 10~ 28 cm 6 s -1 , 



(10) 



2 



where A is due to two-body collisions, and B is due to three-body collisions. The experimental value of the two-body 
loss rate A is consistent with the value give above, provided it is interpreted as the loss rate per event ]lfj| . In view 
of this fact, we reinterpret both A and B as rates per event rather than per particle. Since 2 particles are lost in a 
two-body collision, and 3 particles are lost in a three-body collision, the loss rates per particle is taken to be 2A for 
two-body collisions, and 3B for three-particle collisions UM. 

The gain and loss rates can be taken into account by adding absorptive terms to the NLSE: 
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V 2 + V{r) -t/ H : 
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The initial wave function is taken to be the ground-state eigenfunction in the harmonic potential: 



Mr) = 3/2 0Voexp (-r 2 /d 2 ) 



(11) 
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where No is the initial number of particles. Note that the coefficients A and B were defined as rates rather than 
coupling constants, and do not need corrections factors arising from Bose symmetry. 
We assume that tp is spherically symmetric, and put 

1 

r = -ut, 



do 



N 3/2 U 

4tt p 



The dimensionless NLSE for the reduced wave function u(p,r) reads 



du 



cP_ 
dp 2 
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where 
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\_a\_ 
do 



4.87 x 10" 



7 = - = 6.18 x 10" 



A 



2iruido\a\ 
W 



(4tt) ud^\a\ 



= 1.46 x 10~ 4 , 

2.616 x 10~ 5 . 



The instantaneous particle number, given by 



N(t) 



dp\u(p,r)f 



is no longer a constant of the motion. The initial wave function is chosen to be 

u(p, 0) = 27T- 1 / 4 V^Vopexp (V/2) . 
The initial particle number No ~ N(0) is the only adjustable parameter of the model. 



(15) 



(16) 
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III. COMPARISON WITH EXPERIMENTS 



We solve the modified NLSE (14) numerically, using an algorithm due to Goldberg et al |L^ |, which is a version 
of the popular Crank- Nicholson method E3|. The computational spatial lattice has size Lp = 10000, with spacing 
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dp = 0.0008, which corresponds to 2.53 x 10 -7 cm. The temporal grid is of size L T = 40000, with spacing dr — 0.017, 
which corresponds to 3.74 x 10~ 5 s. The boundary conditions are u(0,r) = u(L p ,t) = 0. 

The Rice experiment Q consists of three groups of measurements of N(t) with different Nq. The initial values 
were quoted as No = 500, 100, 0, with errors of ±60. Each data point was the average of 5 separate runs. We fit the 
data by choosing the best Nq for each group, and also perform an average over a uniform distribution of Nq centered 
about the chosen value, with a width adjusted to yield the best visual fit. The three groups of calculations correspond 
to the following distributions: 

Group (No) Distribution 

I 325 385,355,325,295,265 

II 75 85,80,75,70.65 

III 31 61,46,31.16,1 

The comparison with data is shown in Fig.l. For different Nq, the calculated N(t) has practically the the same 
functional form, except that the time origin is shifted. Typically, N(t) rises exponentially to a maximum of 1200, and 
then suddenly collapses to about 450. Thereafter, it exhibits a periodic saw-tooth pattern of growth and collapse. When 
averaged over a uniform distribution of Nq, the initial rise remains unchanged, but the subsequent cycles tend to be 
smoothed out. The experimental distribution was surely not uniform, but no attempts were made to vary the 
distribution function to achieve a more precise fit. 

The behavior of the initial rise depends mainly on the growth coefficient 7, and is insensitive to the loss coefficients 
a and f3. No discernible difference is found by varying a and (3 by factors of 2 to 3. On the other hand, the behavior 
after the first collapse is sensitive to a, /3, and to variations of the initial wave function. Unfortunately, the data in 
that region is not good enough for comparison with theory. 

In our model, we have neglected quantum tunneling, the saturation of the gain rate, trap asymmetry, and other 
effects. The fact that we can fit the data indicates that these effects are of secondary importance. 



IV. DETAILS OF THE COLLAPSE PROCESS 



We can understand some features of the collapse process by examining the wave function just before and after the 
collapse. In order to mark the time more precisely, we show in Fig. 2 a plot of N as a function of time t, which is 
measured in units of 50 dr = 1.87 x 10~ 3 s. With an initial value No = 325, the first collapse happens at t = 127. 

In Fig. 3, we show the reduced wave functions u(x, t) for t = 100, 127, 128. The position x is the site number x 
of the computational lattice. We see that, as time goes on, the peak of the wave function narrows as it migrates to 
smaller x, reaching x = 300 at t = 127, when the collapse begins. At t = 128 the peak has moved to shown 
in the inset in Fig. 3. The particles under the peak are thus suddenly squeezed into a state of very high density. This 
is in the region referred to as the black hole. The area under the peak, minus background, corresponds to about 750 
partilces. The peak has disapeared by t = 130, as we can see in Fig. 4. This indicates that the loss mechanisms, which 
are proportional to high powers of the density, became activated suddenly, and drained the black hole. Once this 
happens, the density drops, and the loss mechanisms become dormant. The particle number grows due to the gain 
mechanicm, and the wave function slowly recovers an approximately Gaussian shape, as seen in Fig. 4. 

According to Fig.l, 1200 — 450 = 750 particles were lost during the first collapse. This is consistent with loss in 
the balck hole, and indicates that loss mechanisms are effective only in the black hole dring a very brief time interval. 
This explains why the behavior of N(t) is not sensitive to the loss coefficients. In our coarse-grained time, the interval 
of loss amount to 0.004 s. Although this is very short compared to the oscillator period, it is long compared to atomic 
collison times. A study at smaller time scales would require going beyond the mean-field approxiamtion. Work in 
this direction has been initiated by Saito and Ueda |14j . 

As we can see from Fig. 4, after the first collapse, the wave function acquires a persistent chaotic component. This 
means that the wave function it is very sensitive to variations in the parameters of the model, and the initial wave 
function, as has been pointed out In Fig. 5, we examine the chaotic feature moser closely by plotting at different 
times the superfluid velocity field 

v s = (18) 
m or 

where <p is the phase of the condensate wave function. At the onset of collapse at t = 127, the velocity is everywhere 
negative, indicating a flow towards the center. It becomes oscillatory at t = 128, and present a more an more chaotic 
look as time goes on. This suggests that the Li condensate may be a good medium for experimental and theoretical 
studies of superfluid turbulence. 
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Figure Captions 

Fig.l. Grow-collapse cycles of the particle number. Data points are solid circles, with error bars omitted for clarity. 
The statistical error is typically 20%. Thick solid lines are predictions of the model for given No- The dotted lines 
indicate results of averaging over a uniform distribution of Nq described in the text. 

Fig. 2. Particle number as function of. time t measure in units of 50 computational steps, corresponding to 
1.87 x 10~ 3 s. 

Fig. 3. The reduced condensate wave function just before and just after the collapse. The position x is the site 
number on the computational grid, with a spacing corresponding to 2.53 x 10~ 7 cm. The inset shows behaviors near 
the origin, in the black hole. At t = 128, there are about 750 particles under the peak at x = 5. They are purged 
from the system during the next time step, in about 0.002 s. 

Fig. 4. After the collapse, the loss mechanism become inactive, and the system gains particles from the thermal 
cloud. The wave function recovers an approximately Gaussian shape, but a chaotic component remains. 

Fig. 5. The superfluid velocity field at various times. There is a flow towards the center just before the collapse at 
t = 127. Thereafter it looks increasing turbulent. 
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